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ABSTRACT 

In order to understand star formation it is important to understand the dynamics 
of atomic and molecular clouds in the interstellar medium (ISM) . Nonlinear hydrody- 
namic flows are a key component to the ISM. One route by which nonlinear flows arise 
is the onset and evolution of interfacial instabilities. Interfacial instabilities act to mod- 
ify the interface between gas components at different densities and temperatures. Such 
an interface may be subject to a host of instabilities, including the Rayleigh- Taylor, 
Kelvin-Helmholtz, and Richtmyer-Meshkov instabilities. Recently, a new density in- 
terface instability was identified. This self-gravity interfacial instability (SGI) causes 
any displacement of the interface to grow on roughly a free-fall time scale, even when 
the perturbation wavelength is much less than the Jeans length. In previous work, we 
used numerical simulations to confirm the expectations of linear theory and examine 
the nonlinear evolution of the SGI. We now continue our study by generalizing our 
initial conditions to allow the acceleration due to self-gravity to be non-zero across 
the interface. We also consider the behaviour of the SGI for perturbation wavelengths 
near the Jeans wavelength. We conclude that the action of self-gravity across a density 
interface may play a significant role in the ISM either by fueling the growth of new 
instabilities or modifying the evolution of existing instabilities. 
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1 INTRODUCTION 

Hydrodynamic instabilities of interfaces in the interstellar 
medium (ISM) have received renewed attention in part due 
the remarkable Hubble images of "elephant trunks" and "pil- 
lars" in the Eagle Nebula (Hester et al. 1996; Pound 1998). 
Most of the theoretical and simulation work has focused on 
the stability/instability of ionization fronts driven by UV ra- 
diation of nearby OB stars (e.g., Mizuta et al. 2005; Williams 
2002). The fronts may undergo acceleration so as to give 
the classical Rayleigh- Taylor instability but the influence of 
recombination may suppress the instability (Mizuta et al. 
2005). 

Hunter, Whitaker, and Lovelace (1997, 1998; hereafter 
Papers 1 and 2) studied the stability of interfaces in more 
quiescent regions of the ISM where the ionizing radiation 
is not important. They identified a new interfacial insta- 
bility which is driven by self-gravity and acts at a density 
discontinuity. This self-gravity interfacial instability (SGI) 
persists in the static limit for all wavelengths and occurs 
in addition to the classical Rayleigh- Taylor instability. Us- 
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ing a normal mode analysis, they derived the linear growth 
rate of the SGI in compressible media in relative motion (al- 
lowing for the influence of Kelvin-Helmholtz instability). In 
the incompressible limit, the growth rate for a perturbation 
~ exp{—iujt) of a planar interface is 



-2ivG{p2 



Pi) 



P2 + Pi 



+ 



gk{p2 



■Pi) 



p2 + pi 



(1) 



In equation iQ, fc is the horizontal perturbation wave num- 
ber (fe > 0), <; a constant background acceleration, and G 
the gravitational constant. The mass densities of the lower 
and upper fluids are specifled as pi and p2, respectively. The 
second term gives the growth rate for the incompressible 
Rayleigh- Taylor (RT) instability, and the flrst term is the 
incompressible growth rate for the SGI. Both instabilities 
persist in the static limit, but several important differences 
are evident. Self-gravity knows no preferred direction, so the 
SGI is destabilizing across any density interface. An inter- 
face is RT unstable only if g{p2 — Pi) < 0, such that the 
heavy fluid sits "on top" of the light fluid. The SGI growth 
rate depends upon the absolute densities of the fluids and 
their ratio but not the perturbation wavelength. The RT 
growth rate changes with perturbation wavelength and den- 
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Figure 1. Schematic of the equiUbrium configuration. 



sity ratio, but it does not depend upon the absolute densities 
in the fluids. Given these dependencies, the SGI is expected 
to grow faster than the RT instabihty for a fixed value of g 
when the perturbation wavelength is long enough such that 
A = 2-K/k > g / G\p2 — pi\- In planar geometry, the growth 
rate for the SGI depends only weakly upon the Jeans crite- 
rion for the luUy compressible case and not at all upon the 
perturbation wavelength in the incompressible limit. The 
underlying reason that self-gravity is able to drive an insta- 
bility for any wave number is that the configuration is not 
one of minimum energy. 

Numerical studies contrasting the behaviours of the SGI 
and RT instabilities have been performed by Hueckstaedt 
and Hunter (2001; Paper 3) and by Hueckstaedt, Peterson, 
and Hunter (2005; Paper 4). In the nonlinear regime, the 
SGI evolves such that the growth of tenuous bubbles out- 
paces that of dense spikes; whereas, the RT instability is 
characterized by dense spikes streaming into the tenuous 
fluid. In previous work, we sought to isolate the SGI by 
creating a set of hydrostatic initial conditions such that the 
pressure gradient and the acceleration due to self-gravity are 
both zero across the density interface. In the present work, 
we relax this restriction and allow them to be non-zero at 
the interface. The self-gravitational acceleration can effec- 
tively drive an RT-like instability provided the acceleration 
is in the same direction as the density gradient across the 
interface. Although the self-gravitational acceleration varies 
with position, it may be roughly constant near the inter- 
face where the growth is strongest. This RT-like instability 
cannot be strictly isolated because the SGI always will arise 
in the presence of a density interface. Notwithstanding, the 
growth rate of the RT-like instability can be adjusted to 
grow faster than the SGI by increasing the acceleration at 
the interface or decreasing the perturbation wavelength. 

Section 2 describes the envisioned equilibrium config- 
urations, while §3 describes the two-dimensional computer 
simulations. Section 4 discusses the results of our simula- 
tions, and §5 gives the conclusions of this work. 



2 EQUILIBRIUM CONFIGURATION 

In order to isolate the gravitational instabilities, simulations 
are begun from a state of hydrostatic equilibrium. The ge- 
ometry of the problem is shown in Figure The interfacial 



values of the equilibrium densities and temperatures in the 
two regions are pi(0) > P2(0) and Ti(0) < 72(0), and the 
interfacial pressure Pq is the same for both media. The den- 
sities and pressures lapse to zero at heights h\ below and h2 
above the interface. For RT simulations, a constant acceler- 
ation g is applied in the upward direction. 

Adapting a polytropic equation of state (P = Kp^), ex- 
act equilibrium solutions for a self-gravitating gas exist in 
the one dimensional, unperturbed, planar problem when 7 
= 2/3, 1, 2, and 00. A choice of 7 = 2 results in a stiffer 
equation of state than is realistic for the ISM, but it ad- 
mits simple, spatially bounded, analytic expressions for the 
equilibrium distributions. Previous studies revealed no dif- 
ferences in morphology or growth rate betw een simulations 
with 7 = 2 and 7 = 1.4 ijHueckstaedtlbOOll^ . For 7 = 2, the 
hydrostatic and Poisson equations reduce to 



= Knfn(z) 
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dz 
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dz 

In the above, n = 1 for region 1 and n = 2 for region 2, 
Kn = [pn (0)]V2Po, and /„(z) are the self-gravitational ac- 
celerations. The solutions for the densities and accelerations 
have the forms 

and 



P«{z) = Pu{0) cos (^^j + p„(0)^ sin (^-^ 
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The acceleration is continuous across the interface, /i(0) 
/2(0) = fo- The constant /3 is defined as 



/3 = 2V2ttGPo . 

The gravitational scale heights obey the relation 
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Defining hi and h2 as the absolute values of the heights at 
which the density and pressure go to zero in the correspond- 
ing regions, the surface densities are 

JO 

47rG 
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and 



1 — cos 



(I) 



Using equations ||HJ, and Q, it can be verified that 
Gauss' Theorem is satisfied, 

f2{h2) - fi{hi) = -47rG(a2 + ai) . (10) 

By symmetry, /2(/i2) ~ —fi{hi), which leads to the expres- 
sion 

/2(/i2) = -27rG(a2 + ai) . (11) 
Upon integrating equation ^ from z = to 2: = /i2 , we find 
/2(/i2)-/o = -47rGa2 , (12) 
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3 SIMULATIONS 



/o = 47rG'cr2 + /2(/i2) = 47rGa2 - 27rG(a2 + ai) , (13) 



fo = 27rG(cr2 - o-i) . (14) 

Recalling that the densities in media 1 and 2 lapse to zero at 
z — —hi and z — h2, respectively, it follows from equation 
Q that 



tan 



\ h J fo 



(15) 



and 

u..(ii)4. 

Hereafter, we define 9i — hi/h and 02 = h2/l2, both greater 
than zero. The quadrants in which these angles are defined 
depends upon the sign of fo. If fo is greater than zero, 6i is in 
the firs t quadra nt and 02 is in the second quadrant. Defining 
V' — -yZ/J^ + /q , we have the relations sin^i = (3/'tp, cos = 
/o/t/', 02 ~ 11 — 6-1, sin 02 = sin 01, and cos 02 = — cos0i. 
Therefore, 
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Consequently, 

/o = 27rG(cr2 - <Ji) 



47rG/o^ 
47rG^ 



/o 



(17) 
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an identity. The same process can be applied for the case 
/o < 0, with similar results. Therefore, the solutions form a 
consistent set. 

In view of these results, we adopt the following strategy. 
We set the molecular weight to = 2 for both media. We 
specify the densities and temperatures at the interface and 
calculate Po, 13, h, and h- Then, we select /o and compute 
the density distributions from equation ||1J. The pressure 
and temperatures distributions follow from the equation of 
state. The final step is to use equation JSJ to calculate the 
boundary values of /„ (z) for use by the gravity solver. 

Due to discretization of the distributions across the grid, 
a truly static state is not achieved. We deem the setup to 
be sufficiently static if the motions induced in the unper- 
turbed case are negligible compared to any imposed ve- 
locity perturbations. For example, a typical set of initial 
values at the interface is: pi(0) — 10~^''g cm~^, ri(0) = 
20 K, p2(0) = 0.2 X 10-2°gcm-^ and r2(0) = lOOK. 
The a diabatic so und speed for a temperature of 20 K is 
c = a/7Po/pi(0) = 40,800cms"\ We typically use 5% of 
the sound speed (about 2040 cms~^) as the initial pertur- 
bation amplitude. If allowed to run to 5 e— folding times 
(1.06 X lO^'^s), the highest velocities observed throughout 
the grid are less than lOOcms"^. (This represents a conser- 
vative case; most static models show lower velocities.) This 
is more than an order of magnitude lower than the initial 
velocity perturbation amplitude. Deviations from the static 
solutions are not large enough to affect the results of the 
perturbed simulations. 



As a test of the theory and for understanding the 
nonlinear behaviour of the instabilities we have car- 
ried out two-dimensional hydrodynamic simulations using 
CFDLib (Computational Fluid Dynamics Library), which 
was developed at t he Los Alamos National Laboratory 
jKashiwa et alJll994l) . The system of equations we solve is 
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V ■ g = 47rGp 



where Tjk = pSjk + pvjVk is the stress tensor, g is the grav- 
itational acceleration, e is the specific energy of the fluid, 
w — e + p/p is the enthalpy, and the equation of state is 
p ~ {"/ — l)pe. CFDLib is a finite- volume code well suited 
for problems of all flow speeds. The self-gravitational poten- 
tial is solved for in two-dimensions using the MUDPACK 
multigrid code developed at the National Center for Atmo- 
spheric Research (Adams 1989, 1991). Models are run on a 
2D Cartesian mesh of size 257 x 257. Simulations repeated 
on a 513 x 513 grid show a difference in fine scale structure 
but not in growth rate. The normal velocity components of 
the gas are confined by reflective boundary conditions on 
all sides; whereas, the gravitational potential solver uses pe- 
riodic boundary conditions along the side boundaries and 
specified gradient conditions along the top and bottom. 

We induce perturbations along the interface through a 
velocity function of the form 



v{x, z) = vo cos(kx) exp(— fc|z|) 



(21) 



with vo set to 5% of the sound speed in the denser fiuid. The 
perturbation is localized to the interface due to the exponen- 
tial factor in the vertical, ^-direction. We use a velocity per- 
turbation instead of a spatial perturbation for two reasons. 
While our grid resolution is sufficient to determine growth 
rates and obtain a sensible picture of nonlinear structure, it 
is too coarse to impose a spatial perturbation without giving 
rise to spurious instabilities due to the square cell structure. 
Also, without careful consideration, imposing a perturba- 
tion across an interface gives rise to a decaying as well as 
a growing mode. The effect of the decaying mode upon the 
velocity is easily seen and considered in determining growth 
rates. 

For ease of notation, we quote all times in units of the 
e-folding time for the incompressible, linear SGI {g = 0) as 



determined from equation Q, te 



■ 2.11 X 10"s. We 



normalize all growth rates by dividing by the corresponding 
SGI growth rate, losgi = t'^ = 4.73 x lO'^s-^ We define 
two ratios for each model: Q,s for the growth of dense spikes, 
and Q,B for the growth of tenuous bubbles. We can define an 
average normalized growth rate Q,t ~ 0.5(r2s -I- fls)- This 
language is consistent with typical RT descriptions. 

For this study, we define two different classes of mod- 
els. For the first suite of models, we selected a pertur- 
bation wavelength, A — 3.21 x lO^^cm, a value nearly a 
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factor of ten smaller than the Jeans length in the denser 
fluid. We estimate the acceleration that would give rise to 
an RT instability with the same growth rate as the SGI, 
g = ^G(pi(0) - P2(0))A = 2.43 x IQ-^'^cms-^. (The fac- 
tor \/2 is a geometric correction to account for the use of 
a one-dimensional wave vector instead of a two-dimensional 
wave vector. The value of k is reduced by \/2 in going from 
a 2D to a ID wavevector, so the value of g must increase by 
y/2 to arrive at the correct growth rate. By the same logic, 
we also adjust our definition of the Jeans wavelength in §4.3 
to preserve the relationship k\ = 2-k.) We compare the rate 
of perturbation growth for accelerated (/o 7^ 0) interfaces 
against both the pure SGI (/o — 0) and RT (g > 0, no 
self-gravity) cases. Our strategy is to select /o values that 
should lead to growth rates equal to quisGi, where q — 1, 2, 
and 3. We compare the theoretical normalized growth rates 
q to the calculated values Qt- In addition, we look at models 
having /o — — g± 1 x 10~^°cms~^, in order to ascertain the 
critical /o value defining the boundary between stable and 
unstable behaviour. 

We follow with a second set of models designed to in- 
vestigate the SGI as perturbation wavelengths approach the 
Jeans length. We remove the RT component (/o = g — 0) 
and use a value 7 = 1.1 to allow more compression. When 
/o = 0, the initial distributions in density, temperature, and 
acceleration can be determined numerically for any value of 
7 (Hueckstaedt 2001). At question is whether the SGI can 
drive a system which is otherwise marginally Jeans stable 
toward global collapse. 



4 RESULTS 

In Papers 3 and 4, we compared and contrasted the growth 
of pure SGI and RT instabilities. An example of previous re- 
sults is shown in Figure |5| Density contours for SGI and RT 
models with identical growth rates as determined by equa- 
tionQare plotted for times of 2fe and 4te. As is characteristic 
of the SGI, the tenuous fingers grow more rapidly than the 
dense spikes. All of the models presented in this communi- 
cation share the same theoretical e— folding time for a pure 
SGI instability, te = uj~^ = 2.11 x lO^^s. All times are quote 
in units of tg. Figure |3 shows a typical velocity plot used to 
determine numerical growth rates. The logarithmic values of 
the maximum velocity in the spikes and bubbles are plotted 
for the RT instability. Lines are fit to the velocity points 
and slopes determined to represent the linear growth rates 
of both spikes and bubbles. Numerical growth rates for both 
the RT instability and SGI do not exactly match theoreti- 
cal values. Rather, the ratio of computational to theoretical 
growth rate varies with velocity perturbation amplitude and 
density ratio across the interface (Paper 4). 



(a) (b) 
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fcl id) 

Figure 2. Density contours for the SGI [(a) and (c)] and the 
RT instability [(b) and (d)] for t = 2U and t = AU. The grey- 
scales show the density values Xl0~^^. The e— folding time = 
2.11 X lO'^^ s is determined from equation [Tl and is used for all 
figures. 
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Figure 3. Velocity growth for RT instability. The maximum spike 
(x) and bubble (-I-) velocities are plotted versus time using dif- 
ferent values for the normalization velocity. Lines are drawn to 
represent the linear phase growth, with the slopes determining 
the growth rates. 



4.1 Modifying growth rates with non-zero /o 

In order to generalize our study, we now allow non-zero 
values for /q. However, we do so without an imposed con- 
stant acceleration (i.e. no g term). So the RT-like contribu- 
tion to the perturbation growth comes solely from the self- 
gravitational term at the interface. For /□ > 0, an upward 
acceleration drives the interface in a RT-unstable manner. If 
/o < 0, the RT-like component is stabilizing, but since the 



SGI is always destabilizing, perturbations may still grow if 
the RT-like term is relatively small. We investigate the two 
issues with the set of models summarized in Table 1. First, 
we examine the overall growth rates when an unstable RT- 
like term is added to the SGI at the interface. Next we ex- 
amine the marginal case where /o = —g along with a small 
deviation in each direction. 

The first study consists of three models: model fO for 
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Table 1. Model Parameters 



model 


/o 


Q 


Us 


Qb 




Figure # 


fO 


0.0 


1 


0.773 


0.883 


0.828 


2 


f3 


7.277 


2 


1.800 


1.529 


1.665 


6 


f8 


19.41 


3 


3.090 


2.646 


2.868 


9 


deq 


-2.426 











11 


dpi 


-1.426 












dmi 


-3.428 










12 


Units for Table 1: 


/o - 


* 10-10 


cm s" ^ 
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Figure 5. Initial pressure distributions for (a) RT instability with 
g = 2.43x10^1'' cm s-2, (b) model fO, (c) model f3, and (d) model 
f8. The vertical lines indicate the interface locations. 
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Figure 4. Acceleration due to self-gravity (/) for (a) model fO: 
fo = 0,(b) model f3: /o = 7.28 X 10-l°cms-2, and (c) model 18: 
fo = 1.94 X lO-^cms-^. 



fo = 0, the pure SGI; model f3 for /o 
which is three times the value of g that would result in a 
pure RT instability with the same theoretical linear growth 
rate as model fO; and model f8 for fo = 19.4 x lO-^^cms"^, 
with eight times the canonical value of g. The accelerations 
due to self-gravity at time zero for these three models are 
shown in Figure 0] A higher value of fo results in higher 
accelerations throughout the grid. The resultant equilibrium 
pressure distributions for these models (plus the pure RT 
instability) are shown in Figure |S] For the pure SGI, the 
pressure falls in both directions away from the interface. But 
as fo is increased, the acceleration in the less dense medium 
is positive instead of negative. As a result, a monotonically 
increasing pressure is required to form an equilibrium. 

The values of fo for models f3 and f8 are chosen such 
that by equation Q we expect twice and three times the 
growth rate, respectively, of model fO. In Figure |S| we show 
the evolution of model f3. Recall that the times are quoted 
in units of the e— folding time for model fO. A comparison of 
figures|S|and|5|shows that a non-zero value of fo does indeed 
lead to a hybrid sort of structure. The mushroom caps seen 
at later times are broader for model fS than those in the 
pure RT case but slimmer than those in model fO. Also, 
the Kelvin-Helmholtz roU-ups in model fS have assumed a 
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Figure 6. Density contours for the SGI model f3 at (a) 0.9te, (b) 
l.Ste, (c) 2. lie, and (d) 2.7te. 



more relaxed shape than the tight rolls seen in the pure RT. 
Velocity vectors are plotted for model f3 in Figures |7] and |H] 
at different times to illustrate the circulation patterns that 
develop. 

The hybrid structure persists in model f8 (Figure [nj. 
The mushroom caps are nearer in size than those in model 
f3 to the pure RT case, but the Kelvin-Helmholtz roU-ups 
retain their SGTlike relaxed shape. As listed in Table 1, both 
f3 and fS have larger linear growth rates for the dense spikes 
than the bubbles (fis > ^b)- This is shown graphically for 
model f8 in Figure 1101 In this respect, the RT-like growth 
induced by the relatively large values for fo dominate over 
the SGI component during the early phases of growth. For all 
three models, the calculated average growth rate Qt is a bit 
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Figure 7. Velocity map for model f3 at t = 0.9te- Arbitrary 

density contours are overplotted to show the interface position. Figure 9. Density contours for the SGI model f8 at (a) 0.6te, (b) 

l.Ote, (c) 1.4te, and (d) l.SU- 
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Figure 8. Velocity map for model f3 at t = 2. lie 



lower than the predicted q from simplified, incompressible 
theory. The calculated growth rate for model f3 is twice that 
for model fO, as expected. However, the calculated growth 
rate for model f8 is relatively high at nearly 3.5 times that 
of model fO. This larger than expected growth rate may be 
an indication that the amplified RT component grows too 
fast to be affected by the SGI component of the growth. We 
note that the calculated growth rate for spikes in model f8 
{yis = 3.09) is three times the expected spike growth rate 
for a pure RT instability with g — /q. Thus neither the SGI, 
nor deviations from incompressibility, appear to alter the 
rapid RT-like growth when /o is very large. 

We now ask what happens if we apply a negative value 
of fo. In this case, /o(pi — P2) < so the RT-like compo- 




-1.0L 

0.000 0.125 0.250 0,375 0.500 0.625 0.750 0.875 1.000 
time in e-folding times 



Figure 10. Velocity growth for model f8. The spike and bubble 
velocities are plotted with different normalization values. 



nent is stabilizing. The evolution for model deq (Figure lTTll . 
for which fo = —g, is not quite stable even though incom- 
pressible theory predicts 5 = 0. The structure resembles the 
growth of the SGI more than that of the RT instability. The 
inherently unstable SGI seems more robust than the sta- 
bilizing RT component. If the RT component is weakened 
(decreased in absolute value), as in model dpi, the result is 
unchanged except for a slight increase in the rate of growth. 
If the RT component is increased in absolute value (model 
dmi), the SGI growth is quenched (Figure IT^ . 
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Table 2. Model Parameters 




model 


A/Aj 


Pmax (3te ) 


Pmax(6te) 


Figure # 


LJl 


1 


10.8 


20.4 


16 


LJ2 


1/2 


10.7 


12.2 


13 


LJ3 


1/3 


10.5 


11.1 


15 



Aj = 2.92 X 10 cm; densities given in units of 10 gc 
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Figure 11. Density contours for the SGI with /o = —g = —2.43 X 
10-"'cms-2 at (a) 2.0te, (b) 3.0te, (c) i.Ote, and (d) b.Ote- 
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Figure 12. Density contours for the SGI with /o = —3.43 X 
lO-i^cms-^ at (a) 2.0ie, (b) S.Ote, (c) 4.0ie, and (d) 5.0te. 



4.2 Long wavelength perturbations with 7 = 1.1 

For another set of models, we decreased the ratio of spe- 
cific heats and set /q — 0. We seek to examine the be- 
haviour of the SGI as perturbation wavelengths approach 
the Jeans wavelength for gravitational collapse in the rela- 
tively dense, cool gas of medium 1 at the interface. We recall 
that the Jean wavelength (Aj) is the wavelength at which 
the force of self-gravity exactly balances the stabilizing pres- 
sure force in an infinite, uniform, isothermal medium. Ad- 
justing for non-unity 7, the t hree-dimen sional Jeans length 
is given by Aj3d = pT^O) y/ (TvyPo) /G = 2.07 x 10'**cm. 
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Figure 13. Density contours for the SGI model LJ2 with A 
Aj/2 at (a) 3.0te, (b) 4.0te, (c) 5.0te, and (d) 6.0ie. 



We multiply this result by %/2 to account for the two- 
dimensional nature of our simulations and arrive at the value 
A,/ — 2.92 X lO'^^cm. In an idealized model, no gravitational 
collapse is expected when A = Aj; a significantly longer 
wavelength is required for collapse to occur. In contrast, the 
SGI grows for all wavelengths, generating velocities which 
can lead to large compression if given an additional gravita- 
tional boost. 

Other than the change to 7 = 1.1, the physical pa- 
rameters for the three models listed in Table 2 are identical 
to model fO. We increase the cell size in both directions to 
Aj: — Ay = 1.1373 X lO'^^cm, so that the placement of two 
waves across 257 cells results in a perturbation wavelength 
A = Aj/2 (model LJ2). The evolution of this model is shown 
in Figure ITHl The larger cell size results in poorer resolution 
of the roll-up features. This is the price paid to be able to 
examine large scale collapse behaviour without greatly in- 
creasing the computational cost. Velocity vectors for model 
LJ2 are shown in Figure 1141 for t — 5te . By placing three 
waves across the 257 cell extent of the grid, we arrive at 
model LJ3 with A — Aj/3 (Figure As expected, the 
shorter wavelength of model LJ3 results in a growth rate 
and morphology similar to those for model LJ2. 

We observe different behaviour when A = Aj (model 
LJl). For this model, we extended the computational grid 
in the x direction to 513 cells and maintained the same cell 
size. As shown in Figure 1161 the downward motion of the 
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tenuous bubbles (which outpaces the upward moving spikes 
in models LJ2 and LJ3) is slower than observed for the other 
models. Only the central part of the grid is shown in order 
to maintain a visual scale consistent with the previous fig- 
ures. Velocity vectors for model LJl are shown in Figure 
1171 for t = 5te. Although the maximum velocities appear- 
ing in Figures FPfl and [T7I are comparable (about 0.4 km/s), 
model LJl exhibits a greater degree of collapse as indicated 
by the converging velocity vectors. The collapsing nature of 
model LJl is highlighted by taking density line-outs through 
dense columns for all three models at i = 5te (Figure 1181 
and t = 6te (Figure The densest regions in models 

LJ2 and LJ3 are found in the dense mushroom caps which 



move upward over time. In model LJl, the densest features 
form near the interface. At t = 6te, the highest compres- 
sion is seen to have occurred at the stem of the spike rather 
than in the cap. We truncate our analysis at t = 6te due 
to boundary effects from the top and bottom of the compu- 
tational grid. Also, we consider the cylindrical symmetry of 
our two-dimensional calculations unrealistic after the onset 
of (three-dimensional) collapse. As an estimate of the mini- 
mum collapsing mass, we note that the mass of a sphere of 
diameter Aj/2 (the positive phase of the perturbation for 
model LJl) and uniform density pi{0) is about 8.2Mq. 
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Figure 18. Density cuts tlirougli columns at t = b.Ote for three 
different wavelengths. The vertical line indicates the initial loca- 
tion of the interface. 
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Figure 19. Density cuts through columns at t = 6te for three 
different wavelengths. 

5 CONCLUSIONS 

Based upon the work summarized here, we conclude that 
large scale, planar density discontinuities (or near disconti- 
nuities) are inherently unstable on a gravitational collapse 
timescalc in two respects: 1. Any density interface is inher- 
ently unstable to the SGI with a growth rate that is scale 
invariant for incompressible media and nearly scale invari- 
ant compressible media. (In compressible media, some of 



the energy goes into compression at the expense of driving 
fluid displacement.) 2. In general, the acceleration due to 
self-gravity will be non-zero at a density interface, which 
will add a Raylcigh- Taylor component to the instability if 
fo{pi — P2) > 0. Those instabilities cause crcnulations along 
an interface to grow into spikes and bubbles with a growth 
rate and nonlinear structure that depend upon the relative 
strength of the SGI and RT-like components. For wave- 
lengths that approach and exceed the Jeans length in the 
denser medium, structures initiated by intcrfacial instabili- 
ties are likely to undergo continued gravitational collapse. 

An inspection of Hubble Telescope images of interstel- 
lar clouds reveals crcnulatcd interfaces such as the "elephant 
trunks" and "pillars" in the Eagle Nebula (Hester et al. 
1996; Pound 1998). These structures have commonly been 
interpreted in terms of the Rayleigh- Taylor instability of an 
ionization front created by the UV radiation from nearby 
OB stars, but without any consideration of the self-gravity 
of the cloud. However, recent simulations of the ionization 
front dynamics indicate that the Rayleigh- Taylor instabil- 
ity is quenched when hydrogen recombination is included 
(Mizuta et al. 2005) as predicted by Kahn (1958). It is clearly 
of interest to investigate the ionization front stability includ- 
ing self-gravity to see if the SGI instabiliy overcomes the sta- 
bilizing effect of recombination. For this we would need to 
modify the third line of equation (20) to include the different 
energy sources and losses due to the UV absorption, recom- 
ination, and radiative cooling. It is also of interest to extend 
our work to larger computational grids which will allow us 
to follow perturbation growth to longer times. A more am- 
bitious objective is to simulate the SGI in three-dimensions 
at both short and long wavelengths. 

Gravitational processes (both global and intcrfacial) 
clearly have an important role in the evolving ISM. For ex- 
ample, Burkert and Hartmann (2004) have shown that grav- 
itational forces give rise to a variety of structures along the 
edges of finite, self-gravitating sheets in a manner consis- 
tent with observations of local molecular clouds. In the final 
stages of star formation, gravitational forces dominate. But 
long before the final stages of collapse, gravitational forces 
arc important for driving instabilities which convert grav- 
itational energy into fiow kinetic energy, enhance density 
inhomogeneities, and determine the partitioning of energy 
between different length scales. 
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